Non— Equilibrium Dynamics of the Anderson Impurity Model 
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The M-channel Anderson impurity model {A4 = 1,2) is studied in the Kondo limit with a finite 
voltage bias applied to the conduction electron reservoirs. Using the Non-Crossing Approximation 
(NCA), we calculate the local spectral functions, the differential conductance, and susceptibility at 
non-zero bias for symmetric as well as asymmetric coupling of the impurity to the leads. We describe 
an effective procedure to solve the NCA integral equations which enables us to reach temperatures 
far below the Kondo scale. This allows us to study the scaling regime where the conductance depends 
on the bias only via a scaling function. Our results are applicable to both tunnel junctions and to 
point contacts. We present a general formula which allows one to go between the two cases of tunnel 
junctions and point contacts. Comparison is also made between the conformal field theory and the 
NCA conduction electron self-energies in the two channel case. 

PACS numbers: 72.10-d, 72.15.Fk, 72.10.Qm, 63.50.+X 



I. INTRODUCTION 

In recent years, the Kondo model and the Anderson im- 
purity model in its Kondo limit have been investigated 
extensively by use aL numerical renormalizatioijL|flToup 
(NRG) calculationsJilD the|Bethe ansatz method,clo con- 
formal field theory (CFT)fl and auxiliary particle tech- 
niques. In this way a consistent theoretical understand- 
ing of the Kondo effect in equilibrium has emerged. In 
particular, the ground state of the system depends on 
the symmetry group of the conduction electron system: 
If the number of channels, M, is less than the level degen- 
eracy, N, the screening of the local moment at energies 
below the Kondo scale, Tk, leads to a singlet Fermi liq- 
uid ground state with strongly renormalized Fermi liquid 
parameters. If, in contrast, M > N, the ground state 
is M-fold degenerate, leading to a non- vanishing entropy 
at zero temperature T and a characteristic energy depen- 
dence of the density of states, obeying fractional power 
law below the Kondo scale. This is mirrored in an anoma- 
lous (non-Fermi liquid) behavior of the thermodynamics 
as well as transport propertieaJ. 

On the other hand, there has been much less work 
on the non-equilibrium Kondo problem, where the elec- 
tron distribution is not in local equilibrium about the 
Kondo impurity and linear response theory is no longer 
sufficient. Possible new effects in this situation include 
the breaking of time reversal symmetry and the appear- 
ance of a new energy scale like the charge transfer rate 
through a tunneling or point contact. The phenomena of 
tunneling through magnetic impurities has been explored 
since the 1960's when zero bias anomalies (ZBA) ■vpjffi'^ ob- 
served in metal-insulator-metal tunnel junctionsoEl The 
origin of these zero bias anoraalies was understood in 
terms of perturbative theories,! which captured the ba- 
sic phenomena: a logarithmic temperature dependence 
and a Zeeman splitting of the ZBA peak in a finite mag- 
netic field. Although the nriginal theories were quite suc- 
cessful in fitting the datajj they were not able to get to 



what we now know as the low temperature strong cou- 
pling regime of the Kondo problem. In view of theoretical 
advances since that time, it is worthwhile to re-examine 
the non-equilibrium Kondo effect, particularly in the low 
temperature regime. 

In addition, there have been a number of new and 
interesting realizations of the Kondo effect in non- 
equilibrium. With recent advances in sample fabrication, 
it has become possible to see ja. zero bias anomaly caused 
by a single Kondo impurity.E3 Even more intriguing is 
the observation of zero bias anomalies in point contacts 
that exhibit logarithmic temperature dependence at high 
temperatures and power law behavior at low tenmfip.- 
tures but no Zeeman splitting in a magnetic fieldliinlj. 
Such zero bias anomalies may be described by the two- 
channel Kondo model, where the ZBA is caused by elec- 
tron assisted tunneling in two-level systems (TLS), aii 
though other descriptions have been projjosed as well.til 

In the 1980's Zawadowski and VladarEj showed that 
if two-level-systems with sufficiently small energy split- 
tings existed in metals, then one could observe a Kondo 
effect due to the electrons scattering off these TLS. In 
this case the TLS plays the role of a pseudo-spin. One 
state of the two-level-system may be regarded as pseudo- 
spin-up and the other as pseudo-spin-down. An electron 
scattering off the TLS can cause the pseudo-spin state of 
the TLS to change. The electron state also changes, e.g. 
its parity, in the process. This electron-assisted pseudo- 
spin-fiip scattering plays the role of spin-flip scattering 
in the standard case of the magnetic Kondo effect. De- 
tailed analysis, taking into account the different partial 
waves for scattering off the TLS shows that a Kondo 
effect is indeed generated by this electron-assisted tun- 
neling. However, since the true electron spin is conserved 
in scattering from the TLS, there are two kinds or chan- 
nels of electrons. Hence tiie the system will display the 
two-channel Kondo effectllB. Level splitting and multi- 
electron scattering may disrupt the two-channel non- 
Fermi liquid behavior. The stability of the two-channel 
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'or systematic improvements of the 
^ Moreover, the NCA may be gen- 



fixed point against these.perturbations is currently a sub- 
ject of investigationEJOEa. 

In recent years a number of techniques have been 
applied to the non-j£fluilibrium Kondo prohlem: vari- 
ational calcjalationsjlj perturbation theory,L£l equatioH, 
of motion,Eil perturbative functional integmJ methods,E3 
and exactly solvable points of the modeLa One of the 
most powerful techiiij3U£s in this context is the auxil- 
iary boson techniqueOtj. It has two major advantages: 
(i) In its lowest order self-consistertt-|aMproximation, the 
non-crossing approximation (NCA)E3iH3, it yields an ac- 
curate quantitative descriptipH pf the single-channel An- 
derson model in equilibriuiriH3~E2l down to low tempera- 
tures, although it does not capture correctly the Fermi 
liquid regime. The NCA even describes the infrared dy- 
namics of the two-channeL^odel correctly as one ap- 
proaches zero temperature.E3 (ii) The NCA is based on 
a standard self-consistent Feynman propagator expan- 
sion. Therefore, in contrast to exact solution meth- 
ods, it need not rely on special symmetry properties 
which are not always realized in experiments. The for- 
malism also al 
approximation, 
eralized for non-equilibrium cases in a straight-forward 
manner. This has been achi&ved recently for the single- 
channel Anderson modelEa'cZI However, the low tempera- 
ture strong coupling regime of the model was not reached. 

In this article we give a formulation of the NCA away 
from equilibrium which allows for a highly efScient nu- 
merical treatment, so that temperatures well inside the 
low energy scaling regime may be reached. In order to en- 
able other researchers in this field to more readily apply 
this method to related problems we describe the numeri- 
cal implementation of the formalism in some detail. Sub- 
sequently, we study a number of non-equilibrium proper- 
ties of the single- and especially the two-channel Ander- 
son model in the Kondo regime: (1) In linear response we 
study the conductance of the two channel model and the 
one channel model for different spin degeneracies. We use 
both tunnel junction and point contact geometries and 
discuss how to go continuously between the two. These 
results are compared to the bulk resistivity. (2) The non- 
linear response is computed for the same Anderson mod- 
els, and the scaling of the differential conductance at low 
temperatures and voltages is studied. (3) The NCA self- 
energies are compared to those obtained by conformal 
field theory. This sheds light on the question how far the 
equilibrium CFT results for the scaling function are ap- 
plicable to non-equilibrium situations. (4) The effect of 
an asymmetry in the coupling of the impurity to the two 
leads is studied and shown to be consistent with asymme- 
tries of ZBA's observed experimentally. (5) Finally, we 
compute the effect of finite bias on the local (pseudo)spin 
susceptibility. Temperature and voltage scaling is verified 
below Tji, but large differences between the temperature 
and the voltage dependence are found outside of the scal- 
ing regime. 



The paper is organized as follows: In section II the 
one- and two-channel Anderson models and their appli- 
cability to tunnel junctions and point contacts with de- 
fects are discussed. Section III contains the formulation 
of the problem within the NCA and discusses its validity 
for the single- and the two-channel case, respectively. An 
effective method for the solution of the NCA equations 
both for equilibrium and for static non-equilibrium is in- 
troduced. Section IV contains the results for the quanti- 
ties mentioned above, which are discussed in comparison 
with equilibrium CFT solutions and experiments, where 
applicable. All the results are summarized in Section V. 



II. THE SU(N)xSU(M) ANDERSON IMPURITY 
MODEL OUT OF EQUILIBRIUM 

A. The model and physical realizations 

The single-channel (M = 1) and multi-channel (M > 
1) Kondo effects occur when a local iV-fold degenerate 
degree of freedom, a = 1, . . . ,N, is coupled via an ex- 
change interaction to M identical conduction electron 
bands, characterized by a continuous density of states 
and a Fermi surface. For example, the ordinary Kondo 
effect occurs when a magnetic impurity is coupled to con- 
duction electrons via an exchange interaction. The im- 
purity with spin S plays the role of the N — 2S + 1 
degenerate degrees of freedom, and there is only one fla- 
vor or channel of conduction electrons so M = 1. There 
are other physical situations where there are M bands of 
conduction electrons which are not scattering into each 
other. In this case one says that there are M channels 
or flavors of electrons. For the M channel Kondo effect, 
the channel or flavor degree of freedom, r = 1, . . . , M, is 
assumed to be conserved by the exchange coupling. 

Because of the non-canonical commutation relations 
of the spin algebra, this model is not easily accessible 
by standard fleld theoretic methods. Rather than work 
directly with the Kondo model, it is frequently more con- 
venient to work with the corresponding Anderson model. 
Within the Anderson model, each of the possible spin or 
pseudo-spin states, a, is represented by a fermionic parti- 
cle. By convention the operator which creates a ferniion 
in the local level a from a conduction electron in channel 
T is denoted by d^^ . Since each of the d-states created 
by c?J.^, represents a different (pseudo)spin state, only 
one of the states should be occupied at a time. In or- 
der to enforce this constraint, we use the auxiliary boson 
techniqueCJ and write dl^ as d^^ — f^bf, where /o- is a 
fermion operator and br is a boson operator describing 
the unoccupied local d-level. The constraint is then writ- 
ten as the operator identity Q — /l/cr + X^r ^r^r — 1- 
In terms of pseudo-fermion operators /o- and slave boson 
operators bf, the SU(N)xSU(M) Anderson model is 
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k,a,T.a. ^ 

+ UMlb,cl^ + h.c.), (1) 

where fa-, {(J — 1, . . . N) transforms according to SU(N) 
and bf, {f — 1, ■ • ■ M) transforms according to the ad- 
joint representation of SU(M). The first term in Eq. (|I|) 
describes the conduction electron bands with kinetic en- 
ergy ep offset by —eVa due to an apphed voltage. The 
index a is equal to L, R for the left and the right reser- 
voirs, respectively. Note that the two reservoirs do not 
constitute different scattering channels in the sense of 
the multi-channel model, since the reservoir index a is 
not conserved by the Kondo interaction. The second and 
third terms represent the energy of d-states and the hy- 
bridization term, respectively. The constraint term is not 
explicitly written in the Hamiltonian Eq. ([T]). Note that 
the local charge Q commutes with the Hamiltonian. 

As discussed in the introduction, there are a number 
of possible physical realizations of the one-channel non- 
equilibrium Kondo model: magnetic impurities in tunnel 
junctions, tunneling through charge traps, and possibly 
even tunneling through quantum dots. In each of these 
models the d-states introduced in the Anderson model 
have physical meaning. In the case of a transition metal 
magnetic impurities the d-states are literally the atomic 
d-states of the impurity. For a charge trap, the d-states 
are the electronic states for the two possible spin ori- 
entations of the trap. The two-channel model has been 
proposed as a possible scenario for the occurrence of non- 
Fermi liquid behavior in some UeSj^^v fermion compounds 
with cubic crystal symmetry.EilE3u3 In that case the oc- 
cupied d-states correspond to the states of a low-lying 
non-magnetic doublet of the rare earth or actinide atoms, 
while the empty levels, described by bf, constitute an ex- 
cited doublet of local orbitala^. 

On the other hand, for the physical realization in terms 
of two-level-systems, the empty states {b\) do not have 
direct physical meaning. They are introduced as a con- 
struct in representing the pseudo-spins such that the 
channel quantum number r is conserved by the Kondo 
interaction. Via a Schrieffer- Wolff transformatiorO one 
can show that the low energy physics of the Ander- 
son model of Eq. is the same as for the Kondo 
model in the limit when the occupation of the d-states 
rid approaches one. Thus, although we use the Anderson 
model, the results for the low energy physics are expected 
to be the same as for the Kondo model. 



B. The Non— crossing Approximation (NCA) 

1. Validity of the NCA 

In the present context we are interested in the Kondo 
regime of the Anderson model Eq. P), where the low en- 



ergy effective coupling Ja = \Ua\'^/sd between the band 
electrons and the impurity is small, Af{0)Ja <C 1, with 
A/'(0) the band electron density of states per (pseudo-) 
spin and channel. The NCA is a self-consistent conserv- 
ing perturbation expansion for the pseudo-fermion and 
slave boson self-energies to first order in J\f{0)Ja. Con- 
sidering the inverse level degeneracy as an expansion 
parameter, the NCA includes all self-energy diagrams 
up to 0{1/N). The self-energies are then made self- 
consistent by inserting the dressed slave particle prop- 
agators in the rfiDynman diagrams instead of the bare 
propagators.EniSj It is easiliy seen that this amounts to 
the summation of all self-energy diagrams without any 
propagator lines crossing each other, hence the name 
Non-crossing Approximation. 

One may expect that the self-consistent perturbative 
approach is valid as long as the summation of higher or- 
ders in Jq or l/A'^ do not produce additional singulaxitiis 
of perturbation theory. It has recently been showroO'EJ 
that such a singularity does arise in the single-channel 
case (M = 1, N = 2) below the Kondo temperature Tk 
due to the incipient formation of the singlet bound state 
between conduction electrons and the local impurity spin. 
However, around and above Tk and in the Kondo limit 
of the two-channel model {N = 2, M = 2) even downJjp 
the lowest temperatures this singularity is not presents. 
Indeed, the NCA has been very successful in describing 
the single-channel Kondo model except for the appear- 
ance of spurious non-analytic behavior at temperatures 
far below Tk- The spurious low temperature properties 
are due to the fact that the NCA neglects vertex correc- 
tions responsible for restoriiittptiie Fermi liquid behavior 
of the single-channel moddjMo. A qualitatively correct 
description was achievedHj'Ej for the wide temperature 
range from well below Tk (but above the breakdown tem- 
perature of NCA) through the crossover region around 
Tk up to the high temperature regime T > Tk- For the 
multi-channel problem, M > N, the complications of the 
appearance of a spin screened Fermi liquid fixed poipt 
are absent. For this case it has recenthj been shownE2l 
that the NCA does reproduce the exactau low-frequency 
power law behavior of all physical properties involving 
the 4-point slave particle correlation functions, like the 
impurity spectral function and the susceptibilities, 
down zero temperature. Therefore, in the-piulti-channel 
case the NCA is a reliable approximationcS for quantities 
involving Ad (like the non-equilibrium conductance) and 
the susceptibilities, even at the lowest temperatures. 

2- NCA in thermodynamic equilibrium 

The slave boson perturbation expansion is initially for- 
mulated in the grand canonical ensemble, i.e. in the 
enlarged Hilbert space of pseudo-fermion and slave bo- 
son degrees of freedom, with a single chemical potential 
— A for both pseudo-fermions and slave bosons. There- 



3 



fore, standard diagram techniques are valid, including 
Wick's theorem. In a second step, the exact projection 
of the equations-ipato the physical Hilbert space, Q — 1, 
is performedcSEjCH. For a brief review of the projection 
technique we refer the reader to appendix A. 

The equations for the self-energies of the retarded 
Green functions of the pseudo-fermions, G''{uj) = (lu — 
Cd — T,^{uj))~-^, and the slave-bosons, D^{uj) — [lu ~ 
W(uj))^^, constrained to the physical subspace, read 

T,''{uj) = M- I deM{uj - e)/(e - uj)D''{e) (2a) 



W{uj)=N- / deNU^Lo^fU^Lo^CU) , (2b) 

TT 



where V = 7r|J7|W(0), M{uj) = 7V(w)/7V(0) is the bare 
density of states of the band electrons, normalized to its 
value at the Fermi level, and /(w) — 1/(1 + e^'^) is the 
Fermi function. The real and imaginary parts of the self- 
energy are related via Kramers-Kroenig relations, e.g. 



Rel]'^(cj) = ^V J de 



Iml]''(e) 



e — UJ 



(3) 



Taking the imaginary part of Eqs. (^ and defining the 
spectral functions for the slave particles as. 



A{uj) = -ImG"^(tj)/7r = -Iml]''(w) /tt 
B{uj) ^ -lmD'\uj)/Tr = -Imn''(u;) \D''{uj)f/TT 

we arrive at the self-consistent equations 



(4) 



^^'^^ = M- I deM{oj - e)/(e - cu)B{e) (5a) 



^N^l deN{e - u)f{e - u)A{e) . (5b) 

Together with the Kramers-Kroenig relations, Eq. (||), 
Eqs. (H) form a complete set of equations to deter- 
mine the slave particle propagators. However, an ad- 
ditional difficulty arises in the construction of phys- 
ical quantities from the auxiliary particle propaga- 
tors. The local impurity propagator, Gd^ariT — t') = 
-(T'{dd,ar(T)d^ ,^^(t')}), is given by the /- 6 correlation 
function—Thus, its spectral function is calculated within 
NCA asB 



Adiuj) = I y dee-^'^ [A{e + cu)B{e) + A{e)B{e 



(6) 



where 



Z = J dee-f^'lNAie) + MB{e)] 



(7) 



is the canonical partition function of the impurity in the 
physical Hilbert space, Q = 1 (see appendix A). The 



requirement that Z be finite implies that the auxiliary 
particle spectral functipHS vanish exponentially below a 
threshold energy iJo.EZlEa Above the threshold, the spec- 
tral functipas show characteristic power law behavior 
originatingpj from the Anderson orthogonality catastro- 
phe. In Eqs. (H), (0) the Boltzmann factor does not allow 
for a direct numerical evaluation of the integrand at neg- 
ative e if /? = l/ksT is large. It is therefore necessary 
to absorb the Boltzmann factor in the spectral functions 
and find solutions for the functions 



a{u) ^ e-'^'^Aiuj) , b{Lo) = e-^'^Biu) 



(8) 



Using e^'^f{uj) ~ f{—Lu), the equations determining a{uj) 
and b{uj) are easily found from Eqs. (||): 

"^"^^ ^M- f deUiio - e)f{Lo - e)h{e) (9a) 



\G^{^W 



=N^-j deN{e - a-)/(c^ - e)a{e) . (9b) 

The equations for the impurity spectral function and the 
partition function then become 

Ad{uj) ^ ^ J de[Aie + Lu)b{e)+a{e)B{e~Lo)] (10) 



Z = j de [Na{e) 



Mb{e)] 



(11) 



In view of the generalization to non-equilibrium, it is 
instructive to realize that the functions a{uj) and b{(jj) 
are proportional to the Fourier transform of the lesser 
Green functions used in the Keldysh techniqueCJ, 

aiu;) = f G<(c.), G<it-t') = -^{f\t')f{t)) 

b{u) = ^D<{oj), D<{t-t') = i{b\t')b{t)), (12) 

and contain information about the distribution functions 
of the slave particles. Henceforth we will call a{Lo) and 
b{Lo) the 'lesser' functions. Eqs. (||), (^ and (||) form a 
set of self-consistent equations which allow for the con- 
struction of the impurity spectral function Ad- 

A significant simplification of the above procedure can 
be achieved by exploiting that in equilibrium the Eqs. (|) 
and (O) are not independent but linked to each other by 
Eq_(@). Hence, we define new functions A{ijj) and Biuj) 



f{-u;)A{u;) = A{u) , f{-Lo)B{u) = B{lo). 



(13) 



By definition, A{uj) and B{uj) do not have threshold be- 
havior, and the spectral functions as well as the lesser 
functions may easily be extracted from them, i. e. a{uj) = 
/(w)i(cj), B{uj) = f{uj)B{uj). Inserting Eq. (ph into 
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Eqs. (^) one obtains the NCA equations for A(lu) and 

^ 71" J /(-^) 

_Ar£ /deAr(e-^)^^i^;Mz!)A(,). (14b) 



|G'-(^)| 
|G'-(u;)| 



TT 



One can convince oneself that the statistical factors ap- 
pearing in these equations are non-divergent in the zero 
temperature limit for all frequencies u, e. Thus, by solv- 
ing the two Eqs. (|lj) instead of the four Eqs. (||), (||), 
one saves a significant amount of integrations. The equa- 
tions are solved numerically by iteration. After finding 
the solution at an elevated temperature, T is gradually 
decreased. As the starting point of the iterations at any 
given T we take the solution at the respective previous 
temperature value. In appendix A we describe an elegant 
and efhcient implementation of the NCA equations which 
leads to a significant improvement in computational pre- 
cision as well as speed. The proper setup of the discrete 
frequency meshes for the numerical integrations in the 
equilibrium and in the non-equilibrium case is dicussed 
in some detail in appendix B. In this way temperatures 
of 1/1000 Tk and below may be reached without much 
effort. The solutions we obtained fulfill the exact sum 
rules 



nd = N j defie)Ad{e) ^Nj dea{e) = Uf 

J deA,(e) = l-(l-l)rv 

typically to within 0.1% or better, where and Uf 
are the occupation numbers of physical d-particles and 
pseudo-fermions in the impurity level, respectively. 

An important quantity is the self-energy ScC^^) of the 
conduction electrons due to scattering off the Kondo or 
Anderson impurities. In the limit of dilute impurity con- 
centration a; ^ 1, it is proportional to the bulk (linear 
response) resistivity of the system and determines the 
renormalized conduction electron density of states, which 
can be measured in tunneling experiments. Below we will 
calculate Ec(w) within NCA in order to compare with 
the CFT prediction for the resistivity in equilibrium on 
one hand, and to compare the linear response result with 
the zero bias conductance calculated from a generalized 
Landauer-Biittiger formalism (see section III A) on the 
other hand. 

Ec{a;) is defined via the impurity averaged conduction 
electron Green function in momentum space, Gcki^^) — 
[w — Ek— 5]c(<x')]~^. In the dilute limit and for pure s-wave 
scattering, Tic{lo) is momentum independent. 



I]c(w) ~ Xt{uj) 



(15) 



where t{uj) is the local T-matrix for scattering off a single 
impurity. According to the Hamiltonian, Eq. (||), t(u)) is 



given exactly in terms the local d-particle propagator and 
reads, e.g. for scattering across the junction [L ^ R), 

t{Lo) = URUlGd{uo) . (16) 



C. NCA for static non— equilibrium 

If we apply a finite bias V , the system is no longer in 
equilibrium. We cannot expect the simple relation Eq. 
(g) between the lesser and the spectral functions to hold 
in this case. Therefore, the trick with introducing the 
functions A and B cannot be performed. Rather, the 
NCA equations have to be derived by means of 
dard non-equilibrium Green function techniques,! 
and one has to solve the equivalent of Eqs. (|5|) and (||) 
for the non-equilibrium case without any further sim- 
plification. Defining in analogy to the equilibrium case 
^LM = 7r|C/L^ii;pA/'(0), the NCA equations for steady 
state non-equilibrium are 

X ^ [r„7V(w-e + /Xa)/(e-u;-/ia)] 



(17a) 



a=L,R 

\G^{u;)\^ 



N r 

— / ^^^(^) 



[Ta^f{e- UJ - ^J,a)fie- UJ ~ fla)] 



(17b) 



a=L,R 



X ^ [raA/'(cj-e + ^a)/(t^-e + /iQ)] 



a=L,R 
h{L0) 



N 



de a{e) 



X E [rQA/'(e - w - /ia)/(w - e -I- ^q)] . 

a=LM 



(18a) 



(18b) 



If the density of states Af{Lu) were a constant, the only dif- 
ference between the equilibrium and the non-equilibrium 
NCA equations would be the replacement of the Fermi 
function by an effective distribution function F^ff given 
by 

Feffie) = ^/(6 - ^^L) + ^/(e - hr), (19) 

-L tot -L tot 

where Ttot — '^l + ^r- Since our density of states is 
a Gaussian with a width much larger than all the other 
energy scales, \ed\, ^tot^ Tk, this is in fact the only signif- 
icant modification of the NCA equations. Numerically, 
the most crucial modification concerns the integration 
mesh. The proper choice of integration meshes is cen- 
tral to the success of the iteration and is discussed in 
Appendix B. 
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III. CURRENT FORMULAE, CONDUCTANCE 
AND SUSCEPTIBILITIES 

A. Current formulae and conductance 

For the case of tunneling through a Kondo impurity, 
the current is directly related to the impurity Green func- 
tions. In paiiicular, the current in the left or in the right 
lead is givenEJ'L3 by a generalized Landauer-Biittiger for- 
mula, 



(20a) 



lL(y) = -N-Tl J dojN{uj ~ 

X [G<{L0)-Aa{L0)f(u:^yiL)\ 



Ir{V) = N-Tr J duoN{u - ^ir) 

X [G<iij)-A4cj)f{u;-^iR)] 



(20b) 



where is the lesser Green function of the impurity. 
It is obtained from the pseudo-fermion and slave boson 
Green functions via 



1 



G2i^)=Z I dea{e)B{e~uj). 



(21) 



Making use of current conservation, I]^ = Ir, and taking 
the wide band limit, where Af{Lu) is taken to be a con- 
stant, the current may be expressed solely in terms of the 
impurity spectral function 



I{V) = N 



e 2rir, 



dujAd{Lo) 
Mi?)]- 



hrL + Tr 

x[fiLU-flL)^fiLU-tlR)]. (22) 

The NCA is a conserving approximation^ Therefore, 
the currents computed for the left and the right leads 
should be the same when evaluated numerically. We have 
checked the current conservation within NCA and found 
that the two currents agree to within 0.5%, which sets a 
limit to the uncertainty for the average current, I(V) — 
{II + /fl)/2. 

In order to obtain the differential conductance, G{V) = 
dI{V)/dV, we perform the numerical derivative (/(Vi) — 
/(V2))/(Vi - V2), and take it as the value of G{V) 
V ~ {Vi + V2)/2. The numerical error involved in this 
procedure could be reduced to as little as 2%. The zero 
bias conductance (ZBC) is the special case of the above 
equations in the limit of vanishing applied voltage V ^ 0. 
The ZBC for a tunnel junction is thus 



G(0,T) = N 



duj 



dm 

duj 



Adico). (23) 



It will be useful to compare this to the linear response 
bulk resistivity for a small density of impurities in a 
metal. The resisiivity, p, is related to the impurity spec- 
tral function viaEj 



1/(0 = const / du! 



du! 



m, 



(24) 



where the impurity scattering rate is t~^{oj) — 
xULU^Ad{uj). The impurity concentration is denoted by 

X. 

Most of our calculations were done with symmetric 
couplings, Vl = Ffl. However, this is not necessarily the 
case in an experimental situation, especially for tunnel 
junctions. When an Anderson impurity is placed inside 
a tunneling barrier of thickness d, the tunneling matrix 
element Ua depends exponentially on the distance z of 
the impurity from the surface of the the barrier. Also, the 
bare energy level Ed of the impurity will be shifted due 
to the approximately linear in z voltage drop inside the 
barrier. In order to investigate the consequences on the 
non-equilibrium conductance, we also performed evalu- 
ations with asymmetric couplings. For simplcity, and in 
order to keep the total coupling Ttot — Tl + rfl constant, 
we assume a linear dependence of the F^'s on z of the 
form Vl = rtot(l — z/d) , Tr = Ytotz/d. We also mod- 
ify ed according to ed{V) = €d + {V/2){\ - Izjd). The 
latter modification turns out to be insignificant as long 
as V\td\- 



B. Tunnel junctions vs. point contacts 

The above formulae for the currents and conductances 
are valid in a tunnel junction geometry where the current 
must flow through the impurity. In a point contact the 
two leads are joined by a small constriction. A current 
lo will flow through the constriction without the impu- 
rity being present. In fact, the impurity will impede the 
current due to additional scattering in the vicinity of the 
constriction. The question arises whether the effect of 
an impurity in a point contact is the same in magnitude 
but opposite in sign. In Appendix C we derive a gen- 
eral formula for the conductance which allows one to go 
continuously between a clean point contact and a tunnel 
junction. In the limit of a clean point contact, where the 
transmission probabilities are close to unity, we find that 
the change in the conductance due to an impurity in a 
point contact has the same form as for a tunnel junction, 
except for a change in sign. Thus, in clean samples the 
results for the current calculated for the tunnel junction 
apply for point contacts as well, if one subtracts out the 
background current, lo- If lo is ohmic, the conductance 
G{V) is shifted by the constant dIo/dV. Aside from this 
shift and sign difference, the conductance signals of a tun- 
nel junction and a clean point contact will be the same. 



C. Susceptibilities 

The impurity contribution to the dynamic (pseudo-) 
spin suseeptibility is calculated using the standard 
formulaec^Ej from the lesser and the spectral function 
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of the pseudo-fermions. The formula for the imaginary 
part reads 



1 .0 



Imx(w) 



de 



[A(e + uj)a{e)-a{e)A{e-uj)] . (25) 



The real part can be obtained by means of a Kramers- 
Kroenig relation: 



Rex(w) 



(26) 



The static susceptibility Xo = = 0) follows directly 
from this equation. Note that in the two-channel An- 
derson model as possibly realized in TLS's, this suscep- 
tibility is not the magnetic susceptibility. Rather, it is 
probed by a field coupling to the impurity pseudospin, 
e.g. a crystal field breaking the degeneracy of the TLS. 



IV. RESULTS 

A. Conductance for one— and two— channel models 
with symmetric couplings 

Using the formulae discussed in the previous section, 
we now present the results obtained from the numerical 
evaluation of the bulk resistivity and of the conductance 
for symmetric couplings. For the evaluations a Gaus- 
sian conduction electron density of states Af{uj) with half 
width D was used. All calculations were done in the 
Kondo regime for the set of parameters Ed = —0.67D, 
Tl = = 0.151?. In order to make the most direct com- 
parison to experiment, the results for the two-channel 
case have been computed for a point contact, and the 
results for the one channel case have been computed for 
a tunnel junction, except for Fig. ^ where we compare 
the scaling behavior of the nonlinear conductance for the 
one- and two-channel models. 



1. Linear response conductance and resistivity 

The low temperature limit of the linear response con- 
ductance shows power law behavior in temperature. The 
exponent is determined by the symmetry of the underly- 
ing Kondo model. As explained in the discussion of the 
NCA, we expect to get quantitatively correct behavior 
for the two channel model, but not for the one channel 
case. In Fig. |l| we show the zero bias correction to the 
conductance G(0, T) for a two channel Kondo impurity 
(TV = M = 2) in a point contact. The zero bias con- 
ductance does show the expectedQ T^/^ dependence at 
low T. Deviations from this power law start at about 
1/4 Tji, where the Kondo temperature Tr- is determined 
from the data as the width at half maximum of the zero 
bias impurity spectral function. Ad, at the lowest calcu- 
lated T (see Fig. |l|). The slope of the T^/^ behavior 
defines a constant B^: 
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FIG. 1. Temperature dependence of the zero bias conduc- 
tance for the two-channel model (M = 2, A*' = 2) in a point 
contact. The zero bias conductance has T^^'^ dependence for 
T < Tk/4:. This can be used to roughly extract Tk from 
the experimental data. Bs (compare Eq. (27)) is a material 
dependent constant which has been divided out. Therefore, 
the slope of the low T fit (dashed fine) is equal to unity. 



G(o,r)-G(o,o) = Bs^l/^ 



(27) 



which we will use below in interpreting the nonlinear con- 
ductance. 

On the other hand, for the one channel case (M — 
1,N = 2) one expects dependence because of the 
Fermi liquid behavior at low temperatures. As shown 
in Fig. H for a tunnel junction, the NCA as a large N 
expansion is not able to reproduce this power law for 

= 2 at temperatures below Tk- Increasing to A^ = 4 
and A^ = 6, the ZBC develops a hump as a function of 
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FIG. 2. Zero bias conductance for tunneling through a sin- 
gle-channel Anderson impurity (Af = 1, A = 2) vs. temper- 
ature. The conductance for a clean point contact in the pres- 
ence of a single-channel Kondi impurity would be obtained 
by subtraction of this curve from a (constant) background 
conductance. The graph for N — 2 shows an almost linear 
T dependence at low T whereas the curves for spin degener- 
acy A'^ = 4 and A = 6 show non-monotonic behavior. The 
humps are due to the fact that the Kondo peak of the spectral 
function Ad{uj) is shifted away from the Fermi energy ef by 
about Tk- For T > Tk all the curves fall like \og{T /Tk) for 
approximately one decade. 
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temperature. This peak is due to the fact that the Kondo 
resonance is shifted away from the Fermi level for N > 2. 
Although we know of no experimental evidence for such 
humps in zero bias anomalies, similar humps have be 
seen in the magnetic susceptibilities of these systems.Ea 
Note that for = 6, a behavior seems to appear at 
temperatures below the hump. The temperature range 
shown here is above the breakdown temperature of NCA, 
below which a fractional power law G(0, T) — G(O.O) cx 
_j^M/{M+N) ^0^1^ appear. 

For a bulk Kondo system it is impossible to measure 
the zero bias conductance of single impurities. Instead, 
one measures the linear response resistivity, p. In Fig. ^ 
we show the impurity contribution to the resistivity for 
one channel impurities with N = 2,4, 6. Only the N = 6 
curve shows a convex dependence on T. In fact, p seems 
to behave like (1 — const{T /TkY) at-the temperatures 
shown, consistent with a Fermi liquid.t3 For N = 2 there 
is no convex temperature dependence even down to T = 
Q.02Tk- Figures ^ and || also serve to illustrate that the 
zero bias conductance and the bulk resistivity for the 
same kind of Kondo impurities do not necessarily have 
the same temperature dependence. 




0.6 I ' ' ' 1 

0.00 0.05 0.10 0.15 0.20 

T/Tk,N=2 



FIG. 3. Bulk resistivity vs. temperature for the M = 1 
channel model, N = 2,4, 6. Of the three curves only = 6 
has a clear convex shape and falls roughly like at low T. 
The N = 2 graph again shows almost linear T dependence. 
Note that the humps in the conductance for A'' = 4 and = 6 
are not present in the bulk resistivity p. 



2. Nonlinear conductance 

Recently, it has been shown@ that the two-channel 
model exhibits scaling of the nonlinear conductance 
G(y, T) as a function of bias V and T of the fornix 

G{V, T) - G(0, T) = B^Tm{A^) . (28) 

Here, H is a universal scaling function which satisfies 
H{0) = and H{x) cx x'^ for a; > 1, and A, are 
non-universal constants. The exponent rj is 1/2 for the 
two-channel model. This scaling ansatz is motivated by 




(eV/kBT)i/2 




eV/keT 

FIG. 4. Scaling plots of the conductance of point contacts 
in the presence of (a) a two-channel impurity (M = N = 2) 
and (b) a one-channel impurity {M = 1, N = 2). With 
Tl ~ r_R and Be determined from the zero bias conduc- 
tance (compare Eq. (27)), there are no adjustable parameters. 
There are two regimes in these plots. For (eV/kBT)^ < 1.5 
the curves collapse onto a single curve and the rescaled con- 
ductance is proportional to (eV/fcsT)^. For larger (eV/kBT)^ 
the rescaled conductance is linear on these plots. There are 
substantial corrections to scaling even at T small compared to 
Tk ■ At even larger biases this linear behavior rounds off, indi- 
cating the breakdown of scaling. The temperatures are given 
in units of the respective Tk for the two- and the one-channel 
case. 

the scaling of the conduction electron self-energy in the 
variables frequency w and temperature T as obtained-by 
CFT in equilibrium.B Scaling behavior is well knownEJ to 
be present also in the equilibrium properties of the single- 
channel model {M — 1, N — 2). Hence, in the case M = 
1 one may expect a scaling form of the non-equilibrium 
conductance similar to Eq. ( p8| ) as well, however with 
Fermi liquid exponent rj — 2. 

In order to examine whether the scaling ansatz is cor- 
rect, the rescaled conductance is plotted as a function 
of {eV/ksTyK The conductance curves for different T 
should collapse onto a single curve with a linear part for 
not too large and not too small arguments: Very large V 
or T would drive the system out of the scaling regime. A 
collapse indeed occurs for low bias V < T. However, for 
larger bias the slope of the linear part shows T depen- 
dence (see Ref. pj] for more details). This shows that 



8 



there are significant T-dependent corrections to scaling, 
indicating that finite bias V and finite temperature T are 
not equivalent as far as scaling is concerned, although 
both parameters have qualitatively similar effects on the 
conductance. 

Figs. ^ (a) and (b) show the scaling plots for the cases 
M = 2 and M = 1, respectively, with = 2 in both 
cases. Whereas the two-channel case shows the behavior 
described above with the expected exponent i] = 1/2, the 
NCA does not give the correct exponent for the single- 
channel model. In fact, the data show approximate scal- 
ing, however, the exponent rj extracted from the NCA 
data appears to be equal to unity rather than 2. This 
seems to reflect the dominant linear temperature depen- 
dence of the ZBC that the NCA produces in this case. 
This shortcoming is another consequence of the negli- 
gence of singular vertex corrections within the NCA. 

B. Self— energy scaling and comparison to CFT 

Two possible origins for the above mentioned finite- 
T corrections to scaling at low temperatures T are: (i) 
The non-equilibrium state brings about terms in the elec- 
tron self-energy which break scaling of the form given by 
CFT, and (ii) the deviations from scaling in equilibrium 
at finite T have large coefficients, restricting the scaling 
regime to T < Tr-. 

Recpatly, it has been shown for the single-channel 
modeO that in non-equilibrium the conductance indeed 
has terms which explicitly break the scaling behavior. 
Though the coefficients of these terms are small, scaling 
in the ordinary sense is clearly violated, i.e. at temper- 
atures well below a crossover temperature (Tk) not all 
corrections to scaling vanish. It is quite possible that 
the two-channel Kondo model behaves in an analogous 
fashion. 

Within the NCA we here investigate case (ii) by exam- 
ination of the behavior of the self-enjergies in equilibrium. 
As mentioned above, within CFTQ the retarded con- 
duction electron self-energy Sc(w, T) of the two-channel 
particle-hole symmetric Kondo model is found to obey 
scaling of the form 

lmi:,{Lo, T) - ImS,(0, T) = bT^/^H{p^) , (29) 

where H is again a universal scaling function of the kind 
introduced above (77 = 1/2). The constant, 6, is non- 
universal. Within CFT the sign of h the sign depends 
on whether one is on the weak coupling or the strong 
coupling side of the (intermediate coupling) fixed pointB. 
The NCA calculation is on the weak coupling side and 
yields a positive constant h (see below), in agreement 
with CFT. However, a direct quantitative comparison 
with CFT is less straight-forward because the model used 
in the CFT calculation is particle-hole (p-h) symmetric, 
while ours is not: The local level has a finite position 
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FIG. 5. Scaling plot for the imaginary part of the re- 
tarded conduction electron self-energy for a small concentra- 
tion [x = 1%) of the M = 2 channel Anderson impurities in a 
noninteracting metal. Temperatures are given in units of Tk- 
ImSc has a minimum that is shifted to positive frequencies 
due to finite temperature effects. The data are scaled with 
respect to the point (ajmi„, ImE(LJmi,i). For frequencies be- 
low ujmin. the self-energy behaves like {lo/TY^'^ and scales well 
up to frequencies of the order of Tk- However, for positive 
frequencies the self-energy is strongly temperature dependent 
and scaling is less perfect. The parameters of the CFT pre- 
diction for the particle-hole symmetric Kondo model (dashed 
line) have been adjusted so that the slope for negative argu- 
ments matches that of the lowest temperature NCA curve. 

below the Fermi level, while the on-site repulsion is taken 
to be infinite. The scaling plot Fig. ^ shows a very strong 
asymmetry about the point w = 0. The temperature 
dependence is also substantial. The parameters of the 
CFT curve (dashed line) have been adjusted so that the 
slope for the linear part in (w/T)^/^ for negative frequen- 
cies matches the slope of the lowest T NCA curve. The 
asymmetry of the conduction electron self-energy may 
be traced back to the p-h asymmetry of our model. P-h 
asymmetry may very well be present in the experimental 
systemsEj. However, no or only small asymmetries are 
seen in the measurements of the nonlinear conductance. 
This is presumably because in the expression for the non- 
linear current, Eq. (|2^), the frequency is integrated over 
both positive and negative values, —V/2 < u < +V/2, 
thus averaging out the asymmetry. This conjecture is 
supported by our calculation of the conduction electron 
self-energy (Fig. ||), which displays .-strong p-h symme- 
try, and the non-linear conductanceEJ (Fig. |^), which is 
almost p-h symmetric. 

Finally, in the Anderson model we also consider the 
impurity electron self-energy Tidi^^T)- It can easily be 
computed from the d-Green function. The two physical 
self-energies Y,d{ijJ,T) and Sc(ci;,r) are nonlinearly re- 
lated via Eqs. (jlj), ( p^ for a system of dilute impurities 
in an equilibrium situation. It should be noted that the 
nonlinear conductance is directly related to the spectral 
function Ad{ui) via Eq. |2^. We first examine whether 
the impurity self-energy Srf(w, T) shows scaling behavior 
in equilibrium in the variables to and T. The imaginary 
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FIG. 6. Scaling plot for the imaginary part of the impurity 
electron self-energy for the M = 2 channel Anderson model. 
Temperatures are given in units of Tk- For different tempera- 
tures T, Im(Ed(tj,T) -ImEd(0,r))/(cT^/2) is plotted versus 
the square root of the scaled frequency, {lj/TY^'^. The con- 
stant c depends on details of the model. The left parts of 
the curves (u) < 0) obey the anticipated square root behavior 
and scale very well for \uj\ < Tk- For ua/T > the NCA 
curves show a strong T-dependence even for T <^ Tk- This 
is a possible origin of the T dependent slopes of the nonlinear 
conductance curves in Fig. 5. However, for modestly large 
frequencies, e.g. lj/T < 4 the lowest T curves seem to fol- 
low square root behavior, too. The general asymmetry of the 
self-energy is a consequence of the particle-hole asymmetry 
of the Anderson model considered here. 

part of the retarded impurity self-energy is negative 
by causality and its absolute value shows a peak at 
the Fermi level. In Fig. 0, we plot {Im'Ed(u},T) — 
/mErf(0,T))/(cri/2) vs. {uj/TY^^ for different tempera- 
tures T. The constant c is positive and depends on the 
parameters of the model. Obviously, the left parts of the 
curves (lo < 0) obey scaling very well. For arguments of 
{uj/Ty^^ outside the scaling regime uj,T < Tk the curves 
bend, and the self-energies grow roughly logarithmically 
(not shown in the figure). In contrast, the w/T > parts 
of the NCA-curves show a strong T-dependence even for 
T <C Tk, again a consequence of p-h asymmetry. Nev- 
ertheless, at negative frequencies scaling with the correct 
power law is established over a wide range of the scaling 
argument (w/T)^/^, as long as uj,T < Tk, even though 
uj can be much larger than T. It is very plausible that 
this scaling is reflected in the scaling of the nonlinear 
conductance in the arguments bias and temperature asj 
observed in experiments and our numerical evaluation.Ej 
Furthermore, the deviations from scaling at positive fre- 
quencies could be another reason for the observed finite 
T-corrections to scaling of the conductance.c£l 

C. Conductance with asymmetric couplings 

Up to this point we have taken the couplings of the im- 
purity to the conduction bands to be equal, F^ = F^;. As 
mentioned before, especially for a tunnel junction there 



is no reason why this should be the case. The NCA- 
Eqs. (p^), ( |l^ ) are not symmetric in the couplings, that 
is Tl Fjj is not a symmetry of the equations. This 
suggests that the differential conductance signals are not 
symmetric about zero bias if F^ ^ F^;. Indeed, the On- 
sager relations for a two terminal measurement only ap- 
ply to the linear response regime. For nonlinear response 
there is no simple relation between I{V) and I{—V). 
However, interchanging both F^ ^ Tji and V ^ —V 
is a symmetry. It is therefore enough to show only the 
conductances for F^ > 1/2. The curves with F^ < 1/2 
can be obtained from the F^ > 1/2 ones by reflection 
about the y-axis. 

An example of such asymmetric conductance curves 
is shown in Fig. 0. The data is for the two-channel 
model, but the qualitative aspects of asymmetry does 
not depend on the channel number. The constant B^: is 
dependent on the asymmetry, but has been divided out 
for better comparison of the curves. The asymmetry is 
pronounced even for moderate deviations from symmet- 
ric coupling. Asymmetric conductance vs. voltage curves 
similar to those shown in. Fig. ^ have been observed in 
Ta-I-Al tunnel junctions,El where they were plotted as an 
odd in voltage contribution to the differential conduc- 
tance. 
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FIG. 7. Nonlinear conductance for the M = 2 channel 
case for asymmetric coupling, Fz, 7^ F^j = 1 — F^. As ex- 
pected from the asymmetry of the NCA equations (16) and 
(17), the conductance signals show a quite strong asymme- 
try about zero bias even for moderate differences in the cou- 
plings. Asymmetries in the conductance have been observed 
in metal-insulator-metal tunnel junctions. 

D. Dynamic and static susceptibility for the two 
channel model 

Finally, we also show results for the static and dynamic 
susceptibilities with and without finite bias. Although it 
is unlikely that one will be able to measure the suscep- 
tibility for a single impurity, the susceptibility is one of 
the clearest measures of the screening of the impurity by 
electrons. All data shown below is for the two channel 
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FIG. 8. Imaginary part of the dynamic susceptibility 
lmx{i^) (arb. units) for the M — 2 channel Ander- 
son model in equilibrium {V = 0). Temperatures are 
given in units of Tk- For T ^ 0, x{^) behaves as 
Imx(i^) — cisign(aj)[l — C2-\/|<^/Tk|], in agreement with ex- 
act results for the two-channel case. This is in contrast to 
the exact linear behavior in the M — 1 channel model. 

model. The results for the usual Kondo model show dif- 
ferent power laws, but the general effect of the finite bias 
is the same. 

In equilibrium, in the zero temperature limit, the dy- 
namic susceptibility defined in Eq. (|5|) is given by a step 
function of the formtS 



Imx(w) — Cisign(a;) 



I-C2W — 



TV' 



(30) 



The NCA approaches this behavior as the temperature 
is reduced. At finite temperature, the step is broadened, 
as shown in Fig. H, with the extrema located at values 
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FIG. 9. Static susceptibility Xo (arb. units) vs. temper- 
ature at zero and at finite bias V for AI = 2, N — 2. In 
equilibrium, Xo shows the characteristic, expected logarithmic 
divergence as T approaches zero for the two-channel model. 
Out of equilibrium, this divergence is cut off at a temperature 
corresponding to the bias V. The inset shows that Xo falls 
with below this cutoff. For high temperatures T S> Tk, 
Xo falls like 1/T (Curie- Weiss law). 
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FIG. 10. Static susceptibility Xo (arb. units) vs. bias V at 
various temperatures T for M — 2, N = 2. Temperatures are 
given in units of Tk- Xo has a very similar dependence on V 
and T as long as V, T < Tk (scaling regime) . Xo drops like 
for r ~ O.ITk and hke log(V) around Tk- However, for 
large V 3> Tk, Xo falls less rapidly with V than with T, see 
Fig. 11. 

which grow roughly with T^/^. The real part follows by 
a Kramers-Kroenig relation and diverges logarithmically 
for Lo ~* 0, again cut off at finite T . As a consequence, 
the static susceptibility Xo = Rex(a; = 0) diverges loga- 
rithmically as T approaches zero in agreement with . 
Fermi liquid behavior, as has been predicted before.E 
This logarithmic divergence is well reproduced by the 
NCA technique, see Fig. || 

Out of equilibrium, the finite bias serves as another 
low energy cutoff, but in a nontrivial manner. If we look 
at the extrema of the imaginary part of the susceptibility 
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FIG. 11. Product of the static susceptibility Xo and tem- 
perature T (bias V) vs. T {V) on a semi-logarithmic scale for 
M — 2, N = 2. The T dependence shows saturation at high 
temperatures and therefore implies the Curie law, Xo oc 1/T. 
However, the V dependence is linear at large bias, implying 
that Xo falls less rapidly with V than with T, Xo oc log{V)/V . 
The y-axis units are such that a "free pseudo-spin" would cor- 
respond to a constant value of 1/2 (Curie-behavior at large 
temperatures) . 
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at low temperature but finite bias {V > T) we find that 
they are located at smaller absolute values than at the 
corresponding temperature. The logarithmic divergence 
of the Rex(w) is cut off at about V, so that the static 
susceptibility does not diverge logarithmically as T ^ 
anymore. Instead, it approaches a (^-dependent) finite 
value with a quadratic T-dependence (see inset of Fig. 

However, this does not signal the return of Fermi 
liquid behavior for T < V, since we still have y/V/T 
behavior of the conductance for V well below Tk- Fig. ^ 
shows the T-dcpcndcncc of Xo for V = l/lOTfj-. 
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FIG. 12. The impurity spectral function Ad(ijj) for the 
M = 2 channel, N = 2 Anderson model for various biases 
V. The half width at half maximum of the zero bias spec- 
tral function (solid line) is a measure of Tk- As the voltage 
is increased to eV = ksTK the Kondo resonance is reduced 
(dash-dotted line). At very large bias the resonance shows a 
shoulder and eventually splits into two peaks. Increasing the 
temperature would wash out the peak splitting and restore a 
single, much less pronounced peak. 

Similar behavior (i.e. quadratic in V for low V , log- 
arithmic ior T < V < Tk) is observed for the V- 
dependence of the static susceptibility (see Fig. |l^) ; how- 
ever, there is a difference in the dependence on T and V 
in the regime Tk < T,V and T,V < Ttot- For large 
T > Tk at zero bias the static susceptibility behaves like 
l/T, indicating Curie- Weiss behavior. However, for large 
V at low temperature Xo falls less rapidly. The difference 
becomes obvious if we plot Vxo vs. log{V) and Txo vs. 
log(T) as shown in Fig. |Tl|. Whereas the T-dependence 
saturates, indicating the free moment at high tempera- 
tures the y-dependence shows linear behavior, leading to 
Xo ^ blog{V)/V. This stresses again the different con- 
sequences of rising T and V once one has left the scaling 
regime T^V < Tk- The latter becomes most obvious if 
we look at the impurity spectral function Adito) for var- 
ious values of the bias, see Fig. |l^. For zero bias and 
low temperatures there is a sharp resonance with width 
Tk (solid line). Increasing the temperature above Tk 
the peak would broaden at the expense of its height (not 
shown). In contrast, if we keep the temperature low, 
T <^ Tk, and increase the bias, the resonance first de- 



velops a shoulder and then splits into two much broader 
but distinguished peaks (see eV — lOksT and 20kBT). 
Increasing the temperature would eventually wash out 
the peak splitting and restore a single, though much less 
pronounced peak. This difference in behavior at large T 
vs. large V is the reason for the breakdown of scaling of 
the conductance for T or I^ larger than Tk and the dif- 
ferent behavior of the susceptibility at large arguments 
(Fig. 0). 

V. CONCLUSION 

In conclusion, we have described in detail the analyti- 
cal foundations and the numerical implementation of the 
NCA integral equations for the one- and two-channel 
Anderson model out of equilibrium. Our algorithms en- 
abled us to reach lower temperatures than previously ob- 
tained, allowing us to study the physics deep inside the 
scaling regime of the two-channel model. 

In linear response, we computed the conductance for 
tunnel junctions and point contacts as well as the bulk 
resistivity. The two-channel data for both properties 
show T^/"^ behavior in agreement with results obtained 
by other methods. For the single-channel model and 
N — 2 we find dominantly linear behavior below Tk- 
For iV = 6 the bulk resistivity drops with (Fermi 
liquid behavior) at the lowest temperatures considered 
in this work; however, the tunnel junction conductance 
rises with T^ , reaches a maximum below the Kondo tem- 
perature Tk and then falls off logarithmically at higher 
T- This "hump" is associated with the fact that the 
Kondo peak of the impurity spectral function is shifted 
away from the Fermi level for values of > 2. 

If we turn on a finite bias V, the Kondo peak of the 
impurity spectral function first diminishes in height and 
broadens, then splits into two peaks located at the ener- 
gies of the two Fermi levels of the leads at a bias of about 
IOTk- The non-equilibrium conductance is again consis- 
tent with linear behavior in the regime T < V < Tk for 
the single-channel case with M = 1, N = 2- Therefore, 
we can plot the conductance as a function of eV/ksT and 
achieve scaling for modest bias V- Whether similar scal- 
ing of the conductance but with argument (eV/kBT)'^ 
can exist for the case TV = 6 is yet to be determined. 
The tunnel junction conductance falls with for bias 
V < Tk- This is in stark contrast to the hump in 
the T-dependence of the zero bias conductance. If at 
all, scaling seems possible only for temperatures well be- 
low the temperature where the hump occurs. The two 
channel data show scaling with respect to the argument 
{eV/kBTy^^ consistent with conductance measurements 
on clean point contacts. It has to be pointed out, though, 
that the scaling at non-zero bias in the two-channel 
as well as in the single-channel model is only approxi- 
mate. Finite T corrections are observed in the numer- 
ical data (and also in the experimental data) for tem- 
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peratures down to about 1/100 Tr-. This indicates that 
finite bias V and finite temperature T are not equiva- 
lent, although they have qualitatively similar effects on 
the conductance. The scaling of the conduction electron 
self-energy turns out to be worse than that of the non- 
linear conductance. This may be traced back to the lack 
of particle-hole symmetry of our model, which leads to 
asymmetries in the self-energy even at the lowest temper- 
ature. Additionally, there are also strong temperature 
dependent corrections to the square root behavior. 

If we allow for asymmetric couplings to the left and 
right Fermi seas, we observe conductance signals which 
are asymmetric about zero bias. Such features have 
been seen in experiments on metal-insulator-metal tunnel 
junctions. 

Finally, we also calculated the dynamic and static 
(pseudo-) spin susceptibility and discussed the modifica- 
tions due to a finite bias by example of the two channel 
model. The dynamic susceptibility approaches a finite 
step at a; = as T — > 0, leading to a logarithmic diver- 
gence of the static susceptibility in agreement with CFT 
results. A finite bias cuts off this logarithmic divergence. 
In a very similar fashion the temperature cuts off the di- 
vergence as the bias is vanishing. Differences in the bias 
and temperature dependence of the static susceptibility 
appear at high bias and temperature outside of the scal- 
ing regime. 
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APPENDIX A: NUMERICAL 
IMPLEMENTATION OF THE NCA EQUATIONS 



limA^oo (C)Gce' 



,/3A 



lim 



A — >oo 



(A2) 



Note that in this expression the factor Q arising from the 
differentiation d/dC, in the numerator may be dropped for 
any operator O whose expectation value in the subspace 
Q — vanishes (like, e.g., O = d„r{t)d'l^{t') or any other 
physically observable operator on the impurity site). The 
canonical (C) partition function is given by 



Z= lim K^Q)gc(A)] Zq=o 



(A3) 



= Zq^o / dee-f^'iNAie) + MB{e)) , (A4) 



where Zq^q is the partition fimction in the subspace Q — 
0. 

The integrals involved in the NCA equations are 
difficult to. compute because of the singular threshold 
structures of ^(w), B{lo), where the position of the 
threshold energy Eo is a priori not known. In order 
to make the numerical evaluations tractable, we apply 
a time t dependent U(V) gauge transformation simul- 
taneously to the / and b particles according to / ^ 
exp(«Aoi)/, b — *■ exp(iAof)6. This transformation is a 
symmetry of the Anderson model and amounts to a shift 
of the slave particle energy or chemical potential by Ao, 
ljj — + -I- Ao- Note that this shift does not affect any 
physical properties, as seen explicitly, e.g., from Eq. (p^). 
After this energy shift, the spectral and lesser functions 
appearing in the NCA equations read 



aAo(w) 

bxA^) 



Iml]'~(w) 



{Lo-ea + Xo- Rel]'-(tj))2 + (ImE'-(a;))2 
[uj + Xo- ReW{uj))^ + (Imn'-(tj))2 



S<H 

nfM 

(w + Ao - Ren'-(w))2 (Imn'-(w))2 



(A5a) 
(A5b) 

(A6a) 
(A6b) 



Below we briefly review the slave boson projection 
technique and describe an implementation which allows 
for a highly accurate as well as efficient numerical treat- 
ment of the singularities of the spectral functions, which 
arise from the projection. 

The exact projection of the expectation value of any 
operator O onto the physical subspace, Q = l,is achieved 
by first taking the statistical average in the grand canoni- 
cal (GC) ensemble with a chemical potential —A for both 
fermions / and bosons b, and then differentiating w.r.t. 
the fugacity ^ = exp(— /3A) and taking the limit A — s- oo, 



{0)c 



lim — 



(Al) 



In particular, we now have from Eq. (A4) 



dC 



= e-"^" / dee-^%NA{e) + MB{e)) , (A7) 
Zq=o J 

The crucial point about making the numerics efficient 
is that Ao is det erm ined in each iteraiiaa, such that the 
integral in Eq. (A7) is equal to unityEJ'Ej. This defini- 
tion of Ao forces the zero of the auxiliary particle energy 
to coincide with the threshold energy, £'o = in each 
iteration step. Thus, it enables us to define fixed fre- 
quency meshes, which do not change from iteration to 
iteration and at the same time resolve the singular behav- 
ior very well, as described in appendix B. The procedure 
described above leads to a substantial gain in precision 
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and significantly improves the convergence of the itera- 
tions, even though the equation determining Ao must be 
solved durin g ea ch iteration. 

From Eq. (A7) and the definition of the impurity con- 



The a' 



are the limits of integration of 



tribution to the free energy, Fimp(T), exp{~(3Fi 



- imp J 



Z{T)/Zq^o{T), it is seen that Ao determined in the above 
way is just equal to Fimp- This provides a convenient way 
of calculating Fimp (T) directly from the auxiliary particle 
Green functions. 



APPENDIX B: INTEGRATION MESHES FOR 
EQUILIBRIUM AND NON-EQUILIBRIUM NCA 

The various features of the auxiliary particle as well as 
the physical spectral functions are charcterized by energy 
scales, which differ by several orders of magnitude. These 
energy scales are the conduction band width D, the local- 
ized level ed, and the dynamically generated Kondo scale, 
Tk, which is typically of order lO^'^D. Moreover, because 
of the T = 0, V — threshold divergence of the auxiliary 
particle spectral functions, the sharpest features have a 
width given by the temperature, which can be of the or- 
der of IQ-'^D. In non-equilibrium, the bias V appears as 
an additional scale. In the numerical solution of the NCA 
equations, discrete, non-equidistant intregration meshes 
must be set up such that all the features at the various 
energy scales are well resolved. 

These meshes can be generated by mapping the grid 
points Xi of an equidistant mesh onto the non-equidistant 
frequency points LUi by means of an appropriately chosen 
function h{x). In the regions where the very sharp fea- 
tures of the spectral functions and the Fermi function ap- 
pear, i.e. near uj = and ui — ±V/2, respectively, we will 
use a logarithmicly dense mesh. On the other hand, in or- 
der to resolve the relatively broad peak centered around 
the local level e^, the substitution LUi = ed + ctSLTi{xi) will 
be used. 

In general, the entire interval of integration is com- 
posed of L meshes {x-}, i = 1 . . .m, I = 1 . . . L. We map 
these meshes onto the nonuniform frequency meshes {u>l} 
via 



1 , 



(Bl) 



We can now rewrite the integration of an arbitrary func- 
tion k{uj) as an integration over the "equidistant" vari- 
ables {x\}: 



dujk(uj) = / dx 

-oo ^ J ai 



dh\x) 
dx 



k(h{x)) 



4=2 



dx^ 



(B2) 



+ 



Bh^ r)h'- 

{x{)k{h{x\)) + —^{xi^)k{h{xl)) 



dx' 



the different regions of the frequency-axis. To cover the 
whole axis we must have a'+^ = 6'. In equilibrium, 
we can get by with four regions: [~oo, —luj), [—luj,0 — 
ep), [0,w/), [a;/,cxD] , where w/ is an interface frequency 
where two regions of the mesh are matched. {\^d\ — 
V > uji » Tk)- By choosing the functions h'{x') as 
£d+ci tan(a;') in the regions with large absolute frequency 
and as C2exp(a;') in the regions \uj\ < ujj we create large 
mesh point spacings far from e p and exponentially small 
spacings ("logarithmic" mesh) at = 0. Proper adjust- 
ment of constants in the /i' 's is required. The frequency 
mesh point spacing near lu — should be at least 10 times 
smaller than T (and/or V out of equilibrium). Crucial 
for the success of this procedure is the introduction of Ao 
(see Appendix A) in the iteration procedure. Ao shifts the 
peaks of the slave particle functions to the neighborhood 
of a; = in each iteration step. This allows to define a 
fixed frequency grid, which leads to a significant increase 
in computational speed and precision. 

Out of equilibrium the distribution function is a dou- 
ble step function with steps at ztV/2. It turns out that 
in the Kondo limit the slave boson spectral and lesser 
functions show broadened peaks at about the same fre- 
quencies. However, the pseudo-fermion functions behave 
differently. They do not split, but have a single peak 
somewhere between the Fermi level and V/2 that shifts 
not linearly with V. To cope with such behavior we wish 
to have good resolution at zLV/2 and at ep- (The latter 
one is to improve the resolution at the location of the 
peak of the pseudo-fermion functions. Unfortunately, we 
do not know how this location will move with increasing 
V) . To achieve this we let the logarithmic mesh end at 
±V/2 coming from larger/smaller frequencies and choose 
the spacing in between according to the sum of two tanh- 
functions which have their zero shifted to ±y/4, respec- 
tively. We have to choose parameters of these functions, 
so that the mesh spacings at the crucial energies is small 
enough to resolve all features of the integrand. These 
parameters depend on the bias V . They have to be cal- 
culated before the mesh is 'set up' whenever we change 
the potential from one run to the next. However, once 
the mesh is set, we do not have to change it anymore 
during the iterations, because of the same reasons as in 
equilibrium. 

The typical total number of integration points used 
is 200 and 250 for equilibrium and out of equilibrium, 
respectively. Out of equilibrium we need about 50 points 
more for the 'inner' region between ±V/2 at moderate 
bias V < 2QTk- For higher bias we have to introduce 
more points in the inner region. Convergence is achieved 
within 100 - 200 iterations. The CPU time to obtain 
a converged solution on a typical workstation is below 
1 minute for the equilibrium case, and of the order of 
minutes for non-equilibrium. 
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APPENDIX C: GENERAL FORMULA FOR THE 
CONDUCTANCE 

In this appendix we derive Eq. ( p^ for the current 
through a constriction with an impurity. We proceed in 
three stages. First, we introduce our scattering state no- 
tation and review the noninteracting case. Next, we de- 
rive a general formula for scattering from an interacting 
impurity. This is valid for point contacts, tunnel junc- 
tions, and anything in between. Finally, we specialize to 
the case of a clean point contact. 

The geometry we consider consists of perfect left (L) 
and right (i?) leads connected by a central region where 
there is scattering. The scattering states, V'(x), are eigen- 
states of the noninteracting problem. They are labeled 
by their incoming wave vectors, k, where fc^ > corre- 
sponds to a right moving wave and fc^ < corresponds 
to a left moving wave, where z is the direction along the 
length of the leads. For example, a state moving from 
left to right {k^ > 0) has the asymptotic form for z 3> 
of 



V'k(x) 



and for z <C of 



(CI) 



(C2) 



--i\k'^ I -s 



^k'^ (x±) 



k', 



The transverse modes, yk^(xj_), in Eqs. ( |Cl| ) and ( |C2| ) 
are chosen to have unit normalization, and Vz — kz/m is 
the velocity along the length of the leads. It is also un- 
derstood that the energy of the incident and transmitted 
waves are the same, ekj^ + ffc^ — Ck'^ + Cfc^ • 

The current for both the interacting and the noninter- 
acting case may be expressed as a cross-sectional integral 
of the 'lesser' Green function: 



2mi 



5<(x,x';tj) 



(C3) 



For the noninteracting case, this Green function may be 
written in terms of the scattering states as 



dkz 
2tt 



^ 2TrS{co - £;k)Vk(x)V:(x')/k(c^), 



(G4) 



where /k(w) is a Fermi function at chemical potential fj,L 
for kz > and at /i^ for kz < 0. Wc will usually refer to 
these Fermi functions as fL{^) and fji^uj), respectively. 
Using the asymptotic expressions of Eqs. ( |G1[ ) and (C2) 
for the the right moving scattering states and the similar 
ones for the left moving states, Eqs. (|C3|) and (C4) lead 



to the usual Landauer formula for the conductance: 



dEk 
2tt 



1= I ^ y : \t5i\'{fL{Ek)~fBiE,)). (G5) 



k I k 



We now add an impurity which includes an interacting 
term to the Hamiltonian. The coupling of the impurity, 
denoted by 0, to the electrons is given by 



H' = 



dkz 
27r 



E 

ki 



Wkcjck + WCclcQ, 



(G6) 



where k refers to the scattering state of incoming wave 
vector k, not a plane wave state. In Eq. (C6) and in the 
previous equations we have not included spin. The entire 
derivation presented here follows through in the presence 
of a spin (or other) index so long as the self-energy is 
diagonal in that index. This is the case for the Anderson 
model used in this paper. In order to simplify the nota- 
tion, we shall proceed without spin and at the end quote 
the final result when the electron spin is included. 

Using Dyson's equation one can express all of the 
Green functions for the full system, g, in terms of the 
noninteracting Green functions, g°, and the full Green 
function at the impurity, g(0, 0). In particular the Green 
function g<(k, k'), which is used to compute the current, 
is given by 



g<(k,k') =2^<5(fc,-fc;)Jk,kl5<(k) 
+ g°^{k)WC9a{0,0)Wu'g:{k') 
+ 5°(k)VFk*ff<(0,0)W^k'<72(k') 
+ 5°(k)M^k*5.(0,0)I^k'<7<(k')- 



(C7) 



In Eq. (C7) all Green functions have the same energy, 
uj. The self-energy a contains th e m any body interaction 
at the impurity site. Equation (C7) is converted to real 
space using 



9r 



,(x,0) = J ^^^k(x)5.(k,0), 

k 



(C8) 



and the similar relation for the advanced Green function. 
The result for the real space g< is 

5<(x,x') =5r(x,0)(7<5a(0,x') (C9) 

^^{^k"(x)+ff.(x,0)I^k"} 
K 

X5°(k") {V^(x') + I^k*".9a(0,x')}. 

As for the noninteracting case, we wish to evaluate the 
current far into the left and right leads. To do this we 
need the the as ymp totic form of the scattering states 
(Eqs. ( |Cl| ) and (C2)) and the asymptotic form of the re- 
tarded and advanced Green functions, 5j,(£i)(x, 0), which 
we define as 



gr(x,0) 
ffr(0,0) 



-E^k'"^^k.(x,)e+(-)'l'=^l^ 



(CIO) 
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Substituting Eq. (C9) into Eq. (|C3|) for the current, then 
yields 

E /£l^^k"+*k3.(0,0)T^k'f/LM (Cll) 

+ E/ ^<|i^5.(o,o)|V<H 

ki 

^^ = -E /^l*kk"+tk5.(0,0)iyk"|V«H (C12) 
- E /^(l^kk"+ikff.(0,0)M^k'f -'5k,k:()/LH 

-E / |^l<l|t£5.(o,o)|V<H 



Eqs. (IcTTl) and (|cT^) are our most general expressions 
for the current. It is useful to compare them to those for 
the noninteracting case (Eg (|C5|) ). W ithout the terms 
involving i7<, equations (Cll) and ( pl2 ) have exactly the 
same structure as the noninteracting current. The effect 
of the impurity is to change the transmission probabil- 
ity for electrons coming from the left or the right. The 
(T< contains the "scattering-out" of an electron from the 
impurity state. This is a new feature of the interacting 

problem. 

Equations (|Cl|) and ( |cT^ ) are valid for an arbitrary 
scattering potential, including both the tunnel junction 
case and the clean point contact case. We model the 
clean point contact case by a perfect wire. The wire 
will have a conductance equal to e'^/h times the number 
of channels at the Fermi energy. The transmission and 
reflection probabilities for this case are 1 and 0: 



i5ki,k'^ 




,RL 



' k,k" 



,LR 



'k.k" 



(C13) 
(CM) 



In this perfect wire case the scattering states are plane 
waves. The impurity is placed at position, x = a, and 
the overlap matrix elements are 



„ik-a 



A 



(C15) 



where A is the cross-sectional area of the wire. As in 
Eqs. ( pl| ) and (|C^, the L here refers to scattering 
states which start on the left, kz > 0, and R refers to 
those which start on the right, fc^ < 0. The distinction 
between left and right moving is a probably unphysical 
here; however, it is useful to make contact to the tunnel 
junction case. Eq. ( |C15 ) implies that 



'A i\vz\' 



(C16) 



Finally, we define the scattering rate of state k from the 
impurity as 

. \W^\^ fdk .11 \W^\' 



A J 2tt 



2\v,\ A ' 
(C17) 



where the integral is done either over fc^ > or fc^ < for 
A = L, R, respectively. Current conservation requires 
that II — Ir so in computing our final result for the 
current we can take any linear combination of and Iji 
which is convenient: 



\W 



L\2 



\w 



R\2 



= E / ^{h{Lo)-fR{u:)) (C18) 

ki 

^^E fr^^''('^)(/i(^) - /«(^))' 

ki k k 

where Aj^iuj) = — Im(7r(0, 0)/7r is the impurity spectral 
function. 

This is our final result for the number current. The first 
term gives the Sharvin point contact conductance. The 
second term is the correction to the current due to the 
presence of the impuxifai. The expression is the same as 
for a tunnel junction,L3L3 except the sign is reversed. The 
correction to the current once one includes the electron 
spin is 



51 



duj 



E 

k I ,s 



OT^L pi? 



(C19) 



where the only change from Eq. (C18) is that there is a 
sum over the spin, s. If we assume a constant density of 
states and no spin dependence of the matrix elements, we 
obtain the expression Eq. (p^, except for the difference 
in sign between the tunnel junction and point contact 
case. Note that in Eq. (22) the and Tji are defined 



with the density of states divided out. 
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